Bending modes, elastic constants and mechanical stability of graphitic systems 
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The thermodynamic and mechanical properties of graphitic systems are strongly dependent on 
the shear elastic constant C44. Using state-of-the-art density functional calculations, we provide 
the first complete determination of their elastic constants and exfoliation energies. We show that 
stacking misorientations lead to a severe lowering of C44 of at least one order of magnitude. The 
lower exfoliation energy and the lower C44 (more bending modes) suggest that flakes with random 
stacking should be easier to exfoliate than the ones with perfect or rhombohedral stacking. We also 
predict ultralow friction behaviour in turbostratic graphitic systems. 



I. INTRODUCTION 

Graphitic systems are used for a wide variety of in- 
dustrial applications, ranging from lubricant and refrac- 
tory materials to neutron moderators in nuclear fission 
reactors 1 and plasma shields in the next generation of 
fusion reactors^. The recent realization of graphene^ 
(single graphitic layer) and the discovery of its unusual 
electronic properties^ have raised the interest on flake 
graphitic systems as a route to produce graphene samples 
of high quality and in large scald^l 

Despite the technological and scientific importance of 
graphitic systems, the knowledge of their elastic proper- 
ties is unexpectedly poor and new insights are needed. 
The values of the elastic constants describe the mechani- 
cal behaviour^ and are decisive in engineering design to 
avoid material failure. In layered materials, they are even 
more important for the thermodynamic properties due 
to a low- lying branch of acoustic vibrations, the bending 
modes, predicted by Lifshitz^ over fifty years ago. Here 
we show that the shear elastic constant C44 affects the 
mechanism of exfoliation that is relevant for the produc- 
tion of graphene. 

The most reliable experimental studies of the elas- 
tic constants have been obtained by inelastic x-ray 
scattering^ and ultrasonic, sonic resonance, and static 
test methods^. The sample used in the first studjff^l 
was single-crystalline Kish graphite, characterized by an 
extraordinary high degree of ordering, the closest ap- 
proximation to the perfect AB stacking graphite (hex-g) . 
The second study 13 was done using highly oriented py- 
rolitic graphite, the closest approximation to turbostratic 
graphite (turbo-g) where the graphitic layers are ran- 
domly oriented around the c-axis. 

Except for C44 and C13, both studies are in agreement 
within the experimental uncertainties. The C13 value in 
turbo-g was determined only by the less accurate static 
test method and it may be affected by errors. Conversely, 
C44 in turbo-g was determined from the sound velocity 
and its value ranges between 0.18-0.35 GPa 13 , one order 



of magnitude lower than 5.03 GPa found in hex-£pU. 

The discrepancy on C4 4 is at tributed to the existence of 
mobile basal dislocation d 13 * 14 *. After neutron irradiation, 
the elastic constant C44 increases by up to an order of 
magnitude^, suggesting that interstitial defects^ could 
pin the dislocation motions, whence the intrinsic value of 
C44 is measured. A principal difficulty with this expla- 
nation is that interstitial atoms inevitably increase the 
shear resistance between graphitic layers and therefore 
they may increase the C44 value by themselves. 

The aim of this study is to investigate from first prin- 
ciples the elastic constants of graphitic systems with re- 
spect to the stacking misorientations between layers, and 
to describe the key role of the shear elastic constant C44 
on the bending modes (thermal property) and mechanical 
stability. We show that stacking misorientations greatly 
affect C44 and that graphitic systems with perfect (hex- 
g) and random (turbo-g) stacking should be considered 
as two distinct materials described by their own elastic 
and thermodynamic properties. 

This paper is structured as follow. In Sec. [TTJ we give 
a brief summary of the theoretical methods and a dis- 
cussion of the LCAO-S 2 +vdW formalism to include the 
long-range van der Waals (vdW) interactions. The con- 
sequences of shear elastic constant C44 on the bending 



modes and mechanical stability are shown in Sec. Ill A 
and |IIIB[ respectively. 

The elastic constants in the case of high-symmetric sys- 
tems (hexagonal, orthorhombic, rhombohedral and A A 
hexagonal stackings) and for graphitic layers randomly 
oriented around the c-axis (turbostratic stacking) are 
presented in Sec. |III C| and |IIID[ respectively. Finally, 
we summarize and comment on our results (Sec. |IV|). 



II. METHOD 

All the calculations are performed using density func- 
tional theory, within the local density approximation 
scheme (LDA), norm-conserving pseudopotentials^ and 
plane waves with cut-off energy of 150 Ry (abinit 
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package^. The k-point mesh was chosen so that the av- 
erage density corresponds approximately to a 32x32x16 
mesh for hex-g. Energies were converged within 0.05 
meV/atom, and elastic constants within 0.5%. For large 
supercells (turbo-g with more than 50 atoms) we have 
used localized basis-set composed of Gaussian orbitals 
(aimpro code) 18 . The elastic constants calculated by the 
two LDA codes are in agreement within 3% or better. 

The choice of LDA is not fortuitous^ and it was 
dictated by test calculations using the generalized- 
gradient approximation (GGA) within the Perdew- 
Burke-Ernzerhof scheme^. According to GGA, the dis- 
tance between graphitic layers is far too large (4.2 A), 
resulting in a negligible interlayer binding energy and al- 
most vanishing out-of-plane elastic constants (C44, C33). 
For these reasons we have dismissed the use of GGA from 
this studjEHtH] 

Even though LDA yields accurate equilibrium dis- 
tances, due to energetical error compensations, the long- 
range van der Waals (vdW) and more generally the weak 
contribu tion t o the out-of-plane interactions is not well 
describe a 24 * 25 l In order to check the importance of these 
effects on C44 and C33, we have used the LCAO-S 2 +vdW 
formalism to include these specific interactions within 
LDAllD ( usm g fireball code)^. 

This formalism takes into account two major contribu- 
tions. The first one, which we call weak chemical interac- 
tion, is a repulsive energy originating from the overlaps of 
electronic densities between the weakly interacting sub- 
systems. Even though the overlaps are rather small, this 
energy is not negligible. This contribution is evaluated 
proceeding to a second order expansion of the electronic 
wavefunctions with respect to the overlaps. 

The second contribution, which is the vdW itself, orig- 
inating from charge fluctuations, can be seen as the in- 
teraction between electronic dipoles. In the frame of 
the dipolar approximation, we use a second order per- 
turbation theory to describe this contribution. This 
method has originally been tested with success on stan- 
dard graphene-graphene interaction, and more recently 
on a wide range of graphitic materials^. In its current 
stage, the analysis of the internal stresses is not imple- 
mented yet, thus the elastic constants that change the 
in-plane bond lengths (Cn, C12, C13) are overestimated 
by up to 15%. However, for C44 and C33 that describe 
the weak interaction between layers, the internal stresses 
are negligible, and the calculated values are expected to 
be very accurate. 

The elastic constants are determined using two differ- 
ent approaches. The first one uses the response-function, 
implemented in the abinit code, to calculate the second 
derivative of the total energy with respect to the strains. 

The second approach uses the elastic energy density 29 . 
For each elastic constant we have applied 21 strain com- 
ponents Sij to the equilibrium crystal structures (eij were 
typically ranging between 0.01 with increments 0.001). 
The elastic constants are found by fitting the calculated 
energies to a polynomial function in the strains. Both 




FIG. 1: (Color online) Transversal acoustic (bending) mode. 
The bending changes the local stacking between graphitic lay- 
ers. The boxes (a-d) show regions with different slopes and 
stackings. The shear Ad gives the deviations from AB phase 
(perfect stacking). 

approaches yield results in agreement within 0.2%. 

III. RESULTS AND DISCUSSION 

A. Bending modes 

The bending modes are atomic vibrations that can be 
excited even at low temperature and strongly influence 
the thermal properties of layered materials. In 1952, 
Lifshitf^ obtained the following dispersion law for the 
out-of-plane acoustic mode u: 

pxuj 2 (q) = C44 (q 2 x + q 2 y ) + C 33 q 2 z + k (q 2 x + q 2 y ) 2 jc (1) 

where p is the density, c is the interlayer distance, q x ,y,z 
are the wave vectors and n describes the intralayer forces 
characterizing the bending rigidity (1.1 eV)P^. 

The small value of C44 (characteristic of graphitic sys- 
tems, see Table [I]) leads to a predominant contribution of 
the transversal bending modes (q z = 0)) in the phonon 
dispersion curves. These modes are sinusoidal displace- 
ments that propagate along the planes and change the 
local stacking between layers (see Fig. [I]). For nearly flat 
planes the shear stacking Ad is almost zero, whereas, 
for positive (or negative) slope Ad becomes positive (or 
negative). Using trigonometric considerations, it can 
be shown that the maximum value of Ad (see boxes in 
Fig. [I]) is given by: 



A \A + (^) 2 

where a is the amplitude, A is the wavelength and c is 
the interlayer distance. 

The crystal resistance to the stacking shear 29 is pro- 
portional to E oc C44 x Ad 2 , and by lowering C44 more 
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bending modes can be excited at lower temperature. In 
the limit of C44 = 0, graphitic systems approach the 
graphene behaviour, where indeed bending modes (or rip- 
ples) are always presentP^. 

B. Mechanical stability 

By imposing the elastic strain energy as positively 
definite 29 , the stability conditions are given by: 

2^13 < C33 (C11 + C12) Cn, C12, C33, C44 > (3) 

Note that C13 does not affect the stability: a positive (or 
negative) value means that under in-plane compression 
the out-of-plane distance tends to expand (or contract). 

The elastic constants Cn, C12 describe in-plane de- 
formations and they possess the highest values due to 
the strong sp 2 bonding interactions within the graphene 
planes. The elastic constant C33 describes out-of-plane 
compression or expansion and it has always a positive 
value for perfect and stacking misorientations (see Ta- 
ble 0. 

The elastic constant C44 corresponds to a shear be- 
tween graphene layers. Due to the weak interaction be- 
tween planes, the C44 value is the lowest and can be 
positive or negative depending on the stacking misorien- 
tations (see Table [i]). The latter elastic constant is the 
only one that can break the mechanical stability condi- 
tion (i.e. C44 < 0). 

C. Elastic constants in high-symmetric 
graphitic systems 

By imposing a translation vector between graphitic 
layers we found four high-symmetric stackings that cor- 
respond to stationary points on the stacking-fault energy 
surface 31 . To calculate the stacking-fault energy surface 
(for graphitic system is often called corrugation energy 
surface/***"^, we nave used 16- atoms unit cell model of 
eight layers in AB sequence in which the stacking at the 
unit cell boundary is changed by imposing a shear dis- 
placement. This represents an intrinsic stacking fault 
between layers at the unit cell boundaries whereas the 
others remains stacked in the AB sequence. The multi- 
ple layers repeat along the c-axis (8 layers) makes neg- 
ligible self-interaction of the intrinsic fault 31 . The sta- 
tionary points, indicated with square and circle symbols 
in Fig. [2^,, correspond to the four high symmetric struc- 
tures (hexagonal, rhombohedral, orthorhombic and AA 
hexagonal graphite, see Fig. |2|d). The in-plane lattice 
parameter was 2.45 A very close to the experimental 
value 2.463A 12 , with no significant changes among all 
the graphitic structures. 

The origin in Fig. corresponds to hex-g (black 
square symbol), the global minimum of the energy sur- 
face with an interlayer separation of 3.34 A (experimental 
value 3.356 A) 12 . With the exception of the C33 value (29 
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FIG. 2: (a) Stacking-fault energy surface (also called cor- 
rugation energy surface). The square and circle symbols in- 
dicate the stationary points corresponding to the following 
high- symmetric structures; (b) The hexagonal, orthorhombic, 
rhombohedral and AA hexagonal stackings viewed perpendic- 
ular (above), parallel (below) to the c-axis. The energies of 
these structures with respect to hex-g are 1.66 meV/atom for 
ortho-g, 0.10 meV/atom for the rhom-g and 9.29 meV/atom 
for AA hex-g (see text). 



GPa), the calculated elastic constants (given in Table |l|) 
are in agreement within the experimental uncertainties 
found in hex-g 12 . Using the LCAO-S 2 +vdW formalism, 
the out-of-plane elastic constant C33 becomes 42 GPa, 
very close to the experimental value 38.7±7 GPsP^. 

The local minimum at 1 /3(1100) (black circle symbol) 
represents the rhombohedral stacking (rhombo-g). This 
structure possesses the same interlayer separation and 
elastic constants of hex-g (the differences are beyond the 
accuracy of the calculations) with formation energy of 
0.10 meV/atom. This very small energy explains why 
the rhombohedral phase is usually 5-15% intermixed with 
the perfect hexagonal one in natural graphitic flakes. 

The saddle point at 1/6(1100) (grey square symbol) rep- 
resents the orthorhombic stacking (ortho-g). The inter- 
layer separation of the primitive unit cell is 3.37 A with 
formation energy of 1.66 meV/atom. The latter energy 
represents the lowest barrier that has to be overcome 
during the shearing process from an ideal configuration 
to another equivalent one. The C44 values are found to 
range between -2.7 GPa, for shearing along (1100), and 
7.7 GPa, for shearing along the (2110) axis. No signif- 
icant changes are found for the other elastic constants 
(see Table [I]). Although unstable (due to the negative 
C44), this structure acts as an intermediate phase during 
the transformation from graphite to diamond 35 . 

The global maximum at 2 /3(1100) (grey circle symbol) 
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FIG. 3: (Color online) (a,b) Representation of accidental corn- 
mensuration. A supercell of vector (n, m) (blue colour) be- 
comes commensurate when rotated by an angle with respect 
to the starting supercell (as we increase n, m the supercell 
surface and number of atoms rapidly increase), (c) Stack- 
ing fault energy surface (or corrugation energy surface) of a 
graphene bilayer with supercell vector (2,1). The respective 
elastic constants are calculated on the local minima of the 
energy surface. Notice that the energy surface becomes much 
flatter leading to a reduction of Caa- 



corresponds to AA hexagonal stacking (AA hex-g). This 
phase has the largest interlayer separation (3.60 A) and 
the highest formation energy (9.29 meV/atom). Its elas- 
tic constants are smaller than in hex-g, and C44 in partic- 
ular, becomes negative (-3.8 GPa) breaking the stability 
conditions. Even though this structure is highly unsta- 
ble, a recent study has suggested that screw dislocations 
locally encourage this stacking 36 . 

In the following section we describe the elastic con- 
stants for graphitic layers randomly oriented around the 
c-axis. 



D. Elastic constants in turbostratic stacking 

The modelling of turbostratic staking is challenging 
since the incommensurate nature of these stackings must 
combine with the finite-size constraint required by calcu- 
lations. To overcome this difficulty we used the method 
proposed by Kolmogorov and CrespP^ in which a layer 
with supercell basis vector (n, m) becomes commensurate 
with a second layer for a specific rotation angle of: 



$ — cos 1 



2n 2 + 2nm — m 2 
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with n > m (4) 



Figure [3^i,b shows the case of the two smaller super- 
cell with basis vector (2,1) and (3, 1) (corresponding to 
14, 26 atoms/layer and rotation angles of 38.21°, 32.20°, 



FIG. 4: An example of 5-layers turbostratic stacking (super- 
cell basis vector 2,1). Each layer is stacked along the c-axis, 
rotated with respect to each other with an angle of 38.21° and 
randomly translated along the basal plane. 



respectively) . Figure [3^ shows the stacking energy sur- 
face for a bilayer supercell of basis vector (2, 1). To bet- 
ter describe the complete misorientation of turbostratic 
stacking we have increased the number of layers along 
the c-axis. Each layer was rotated with respect to each 
other and randomly translated along the basal plane (see 
Fig. [4]). The rotation angles used are 15 values ranging 
from 6.01° to 53.99°. The smallest supercell contains 28 
atoms with 2 layers rotated with respect to each other by 
38.21°, whereas the largest one contains 456 atoms with 
12 layers and rotation angle of 46.83°. 

For each model we carried out extensive structural 
optimizations starting from different translation vec- 
tors along the basal plane. The corrugation energy is 
about one order of magnitude lower than commensu- 
rate structure (see Fig.[2k) with a maximum value of 
0.9 meV/atom (see Fig. [3p). As we increased the size of 
the supercell over the basal plane we found that the av- 
erage corrugation energy tends to decrease up to 20% for 
the largest size (basis vector 8,3 with 194 atoms/layer). 
Extrapolating these results in the ideal case of infinite 
layers, we suggest corrugation energy virtually flat with 
layers mutual independent. 

For all the supercell studied the in-plane lattice param- 
eters remain almost equal to the value found in hex-g, 
whereas the interlayer distances were on average slightly 
larger, 3.42±0.01 A, with formations energies of 3.03=b0.6 
meV/atom. 

The interlayer binding energy between graphitic lay- 
ers (i.e. exfoliation energy) was 21±1 meV/atom (70±4 
meV/atom with vdW), a value slightly lower than the 
24 meV/atom found in hex-g (80 meV/atom with vdW). 
Note that LDA values yield to a binding energy within a 
factor of 2-3 with respect to LCAO-S 2 +vdW formalism 
and experiment values (43 meV/atom found in heat-of- 
wetting experiment!^, 35±10 meV/atom found by an- 
alyzing TEM images of twisted collapsed nanotubes^, 
and 52±5 meV/atom by studying thermal desorption of 
polyaromatic hydrocarbons 39 ). In Table [i] we report the 
calculated values of the elastic constants. With the ex- 
ception of C44, these values hardly change among all the 
studied turbostratic stackings with no clear dependence 
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FIG. 5: (Color online) In-plane lattice parameters vs. pres- 
sure. The solid line represents the results found for hex-g and 
rhombo-g (perfect and rhombohedral stackings). The dashed 
line shows the results found here for turbostratic graphite. 
For comparison, the experimental results are also plotted. 



on the rotation angles and number of layers. The elastic 
constants Cn, C12 slightly decrease by about 3% with re- 
spect to hex-g, remaining within the experimental uncer- 
tainties found in turbostratic samples^. As previously 
found in hex-g, LDA calculations underestimate the C33 
value (27±2 GPa) with respect to the vdW correction 
(36±1 GPa) and the experimental value of 36.5±lGPgP. 

Conversely, we found Ci3=-2.7±0.5 GPa in disagree- 
ment with the experimental value 15±5 GPgP. The lat- 
ter value was only indirectly obtained as a function of the 
other elastic constants by the less accurate static test 
method. This method requires larger strains than ul- 
trasonic experiments and non-linear behaviour of stress- 
strain curve may affect the measured value. Furthermore, 
the linear bulk modulus B a calculated from these elastic 
constants is far too large (2080 GPa), almost double than 
the one found in diamond (1326 GPa^P. 

The linear bulk modulus B a describes the variation of 
the lattice parameter a as a function of the hydrostatic 
pressure 4 - 



41 and it is given by: 



B a 



C33 (Cn + C12) 



C33 — c 



13 



(5) 



This modulus is strongly weighted by C13. For example, 
if we use the measured value found in hex-g 0±3 GPa 
(instead of 15±5 GPa), B a becomes 1240 GPa (instead 
of 2080 GPa). 

Several x-ray studies have measured the linear bulk 
modulus B a . In these experiments powder samples were 
prepared by grinding different types of graphitic materi- 
als, ranging from well crystallized to poorly crystallized 
grain d 42 " 44 [ The good agreement between the experi- 
ments (see Fig. [5| indicates that the linear bulk modu- 
lus B a does not strongly depend on the stacking order 
with a measured value of 1250±70 GPsf^. By includ- 
ing the latter modulus in the Eq. [5j the elastic constant 
Cis becomes 0.3 GPa. Therefore, we conclude that the 
Cis value does not significantly change between turbo-g 
and hex-g and we propose that the same value 0±3 GPa 
should be appropriate also for turbostratic stacking. 
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FIG. 6: (Color online) The elastic constants C44 and C33 as 
a function of the commensurate area over the basal plane in 
turbostratic stacking. The C44 value tends to decrease with 
respect to the supercell size, whereas C33 is not quantitatively 
affected. 



The C44 represents the second derivative of the total 
energy as a function of shear displacement over the basal 
plane. As the corrugation energy of perfect AB stack- 
ing is much higher than turbostratic we are expecting a 
corresponding lower value of the shear elastic constant. 
We have found that C44 tends to decrease as a function 
of the supercell area over the basal plane (see Fig. [6| 
and is independent on the relative rotation angles and 
number of layers along the c-axis (the respective average 
corrugation energies and therefore the C44 values are also 
independent on the relative rotation angles and number 
of layers). 

For commensurate structures of finite area ranging be- 
tween 36-563 A the respective C44 values are within 
0.16-0.31 GPa, one order of magnitude lower than hex-g 
(4.5 GPa) and close to the experimental measures rang- 
ing from 0.18 to 0.35 GPa 13 . The vdW correction does 
not significantly affect the C44 values with respect to 
LDA (see Table l|, suggesting that the variations in en- 
ergy under interlayer shears are nearly identical between 
these two approximations. As we increase the sizes of 
the commensurate structures towards the ideal case of 
infinite layers (incommensurate) the C44 values tend to 
zero with corrugation energy virtually fiat. These results 
indicate that turbostratic stacking possesses the lowest 
friction among all the graphitic materials bearing great 
potential applications in nano-mechanical systems^ 



IV. CONCLUSIONS 

To summarize, we have discussed the importance of 
C44 as the main parameter that restrains the bending 
modes and controls the mechanical stability of layered 
materials. Using advanced ab-initio method, which in- 
cludes vdW interactions, we have provided the first com- 
plete description of the elastic constants in graphitic sys- 
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TABLE I: Elastic constants in unit of GPa for the different graphitic systems. The values between brackets are calculated 
using the LCAO-S 2 +vdW formalism. These results show that the C13 values do not significantly change between turbo-g and 
hex-g and we propose that the same value 0±3 GPa should be appropriate also for turbostratic stacking. The shear elastic 
constants C44 found in turbostratic stacking correspond to commensurate structures of area ranging between 36-563 A . 
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4.4 (4.8) 


-2.7/7.7 (-2.9/7.3) 


-3.8 (-3.8) 



terns. The lower exfoliation energy (3-8 meV/atom) and 
the lower C44 (at least one order of magnitude) found in 
turbostratic stacking suggest that the exfoliation mech- 
anism, relevant for the production of graphene, should 
be easier for graphite flakes with random stacking. Our 
results indicate that turbostratic graphitic systems pos- 
sess the lowest friction among all the graphitic stackings. 
It would be interesting to check these predictions exper- 
imentally. 
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